####### Leslie Finger and Adela Soliz
####### Code From "But Everybody's Doing It!" The Power of Peer ####### Effects in High Peer Density Areas


##### Libraries ############################################################
############################################################################
library(tcltk)
library(MatchIt)
library(cem)
library(Zelig)
library(foreign)
library(AER)
library(xtable)
library(car)
library(mgcv)
library(Matrix)
library(lme4)
library(plm)

##### Function for clustered standard errors ############################################################
se.func <- function(flm,data){
	mat <- estfun(flm)
	dim(mat)
	N <- nrow(mat)
	u <- apply(mat, 2, function(x) tapply(x, data$village, sum))
	u <- na.omit(u)
	vcovCL <- sandwich(flm, meat=crossprod(u)/N)
	ses <- coeftest(flm, vcovCL)
	out <- return(ses)
}


##### Loading data #########################################################
dat <- read.dta("pset2_data.dta")
head(dat)





##########
####Creating individual ids for children
datid <- transform(dat, id=paste(folio, numero, sep=''))
head(dat)
dat <- cbind(datid$id, datid[,3:31])
ncol(dat)
colnames(dat) <- c("id", colnames(dat[2:30]))
head(dat)
ncol(dat)


####FUNCTIONS FOR MANIPULATING DATA

### Getting rid of duplicates to figure out fraction eligible both years.  Then fraction eligible in 98
datdup<- subset(dat, select=c(id,eligible,village))
datdup2<- subset(datdup, !duplicated(datdup))
head(datdup2) 
nrow(datdup2)
fracfunction <- function(x){
	s1 <- sum(datdup2[datdup2$village == x,]$eligible == "pobre")
	s2 <- sum(datdup2[datdup2$village == x,]$eligible == "no pobre")
	return(s1/(s1+s2))
length(unique(datdup2$village))}
fracfunction98 <- function(x){
	s1 <- sum(dat[dat$village == x & dat$year==98,]$eligible == "pobre")
	s2 <- sum(dat[dat$village == x & dat$year==98,]$eligible == "no pobre")
	return(s1/(s1+s2))}



###Subsetting to adding proportion ineligible in school
funcinelig97<-function(x){mean(dat97[dat97$village == x & dat97$eligible=="no pobre",]$scbase)}
funcinelig98 <- function(x){mean(dat[dat$village == x & dat$year==98 & dat$eligible=="no pobre",]$sc)}
funcinelig99 <- function(x){mean(dat[dat$village == x & dat$year==99 & dat$eligible=="no pobre",]$sc)}


	
####Getting rid of duplicates to add total proportion enrolled in 97.  Then functions for proportion in school by year
dat97 <- subset(dat, select=c(id, village, scbase, eligible))
head(dat97)
dat97 <- subset(dat97, !duplicated(dat97))
funcsc97 <- function(x){mean(dat97[dat97$village == x,]$scbase)}
funcsc98 <- function(x){mean(dat[dat$village == x & dat$year==98,]$sc)}
funcsc99 <- function(x){mean(dat[dat$village==x & dat$year==99,]$sc)}


####Adding "size"
sizefunction <- function(x){length(dat$village[dat$village==x])}


###APPLYING THE FUNCTIONS
####Fraction eligible over both years (no repeats) and in just 98
frac <- sapply(unique(datdup2$village),fracfunction)
frac98 <- sapply(unique(datdup2$village),fracfunction98)

####Percent in school by village each year
sc97 <- sapply(unique(dat97$village), funcsc97)
sc98 <- sapply(unique(dat$village),funcsc98)
sc99 <-sapply(unique(dat$village),funcsc99)

####Percent ineligible in school by year
inelig97<- sapply(unique(dat97$village), funcinelig97)
inelig98 <- sapply(unique(dat$village),funcinelig98)
inelig99 <-sapply(unique(dat$village),funcinelig99)

####Size of village
size <- sapply(unique(dat$village),sizefunction)



######PUTTING IT ALL TOGETHER FOR VILLAGE DATA SET ##############
head(datvil)
datvil1 <- data.frame(unique(dat$village),frac, frac98, sc97, sc98, sc99, inelig97, inelig98, inelig99, size)
colnames(datvil1) <- c("village",'frac', 'frac98',"sc97", 'sc98','sc99',"inelig97","inelig98","inelig99", 'size')
datvil1$inelig97[datvil1$inelig97=="NaN"]<-0
datvil1$inelig98[datvil1$inelig98=="NaN"]<-0
datvil1$inelig99[datvil1$inelig99=="NaN"]<-0
datTrpool <- subset(dat, select=c(village, Trpool))
datTrpool <- subset(datTrpool, !duplicated(datTrpool))
datvil <- merge(datvil1,datTrpool, by.x="village", by.y="village", all.x=TRUE)


head(datvil)


####MERGING IT FOR INDIVIDUAL DATA SET
dat1 <- merge(dat,datvil1,by.x='village',by.y='village')
head(dat1)
summary(dat1)
nrow(dat1)

####creating vector with state indicators
dat1$state <- rep(0,nrow(dat1))
dat1$state[dat1$hidalgo==0 & dat1$mich==0 & dat1$puebla==0 & dat1$sanluis==0 & dat1$veracruz==0]<-1 #### 1 is Guerrero
dat1$state[dat1$hidalgo==1]<-2
dat1$state[dat1$mich==1]<-3
dat1$state[dat1$puebla==1]<-4
dat1$state[dat1$sanluis==1]<-5
dat1$state[dat1$veracruz==1]<-6
unique(dat1$state)
head(dat1)

head(datvil)

####Columns for different fractions eligible
###fraction eligible by .1
dat1$fraction1 <- rep(0,nrow(dat1))
dat1$fraction1[dat1$frac>=0 & dat1$frac<=.1]<-1
dat1$fraction2 <- rep(0,nrow(dat1))
dat1$fraction2[dat1$frac>=.1 & dat1$frac<=.2]<-1
dat1$fraction3 <- rep(0,nrow(dat1))
dat1$fraction3[dat1$frac>=.2 & dat1$frac<=.3]<-1
dat1$fraction4 <- rep(0,nrow(dat1))
dat1$fraction4[dat1$frac>=.3 & dat1$frac<=.4]<-1
dat1$fraction5 <- rep(0,nrow(dat1))
dat1$fraction5[dat1$frac>=.4 & dat1$frac<=.5]<-1
dat1$fraction6 <- rep(0,nrow(dat1))
dat1$fraction6[dat1$frac>=.5 & dat1$frac<=.6]<-1
dat1$fraction7 <- rep(0,nrow(dat1))
dat1$fraction7[dat1$frac>=.6 & dat1$frac<=.7]<-1
dat1$fraction8 <- rep(0,nrow(dat1))
dat1$fraction8[dat1$frac>=.7 & dat1$frac<=.8]<-1
dat1$fraction9 <- rep(0,nrow(dat1))
dat1$fraction9[dat1$frac>=.8 & dat1$frac<=.9]<-1
dat1$fraction10 <- rep(0,nrow(dat1))
dat1$fraction10[dat1$frac>=.9 & dat1$frac<=1]<-1



###fraction eligible by .2
dat1$fraction1.2 <- rep(0,nrow(dat1))
dat1$fraction1.2[dat1$frac>=0 & dat1$frac<=.2]<-1

dat1$fraction2.2 <- rep(0,nrow(dat1))
dat1$fraction2.2[dat1$frac>=.2 & dat1$frac<=.4]<-1

dat1$fraction3.2 <- rep(0,nrow(dat1))
dat1$fraction3.2[dat1$frac>=.4 & dat1$frac<=.6]<-1

dat1$fraction4.2 <- rep(0,nrow(dat1))
dat1$fraction4.2[dat1$frac>=.6 & dat1$frac<=.8]<-1

dat1$fraction5.2 <- rep(0,nrow(dat1))
dat1$fraction5.2[dat1$frac>=.8 & dat1$frac<=1]<-1
head(dat1)










##### Subsetting for ineligible kids
dat1ineligible<-dat1[dat1$eligible=="no pobre",]

head(dat1ineligible)

## Sub-setting below and above the .67 of fraction eligible
lowelig.5<-dat1ineligible[dat1ineligible$frac<=.6667,]
highelig.5 <-dat1ineligible[dat1ineligible$frac>=.6667,]
summary(lowelig.5$frac)
summary(highelig.5$frac)
nrow(lowelig.5)  ###2153
nrow(highelig.5) ###2058



###Creating interaction terms
dat1ineligible$interact1 <- dat1ineligible$Trpool*dat1ineligible$fraction1.2
dat1ineligible$interact2 <- dat1ineligible$Trpool*dat1ineligible$fraction2.2
dat1ineligible$interact3 <- dat1ineligible$Trpool*dat1ineligible$fraction3.2
dat1ineligible$interact4 <- dat1ineligible$Trpool*dat1ineligible$fraction4.2
dat1ineligible$interact5 <- dat1ineligible$Trpool*dat1ineligible$fraction5.2


###W/ random effects
lm.re <- lm(sc~Trpool+interact1+interact2+interact3+interact4+interact5+gender+indig+welfare+hohedu+hohage+ hohsex+fam_n++dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+fraction1.2+fraction2.2+fraction3.2+fraction4.2, data=dat1ineligible)
summary(lm.re)
se.func(lm.re, dat1ineligible)


###W/ random effects (state and village fixed)
lm.re2 <- lm(sc~Trpool+interact1+interact2+interact3+interact4+interact5+gender+indig+welfare+hohedu+hohage+ hohsex+fam_n++dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+fraction1.2+fraction2.2+fraction3.2+fraction4.2+factor(state)-1, data=dat1ineligible)
se.func(lm.re2, dat1ineligible)


interact <- dat1ineligible$Trpool*dat1ineligible$frac

###W/a continuous interaction
lm2 <- lm(sc~Trpool+interact+gender+indig+welfare+hohedu+hohage+ hohsex+fam_n++dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r, data=dat1ineligible)
summary(lm2)
se.func(lm2, dat1ineligible)







## IV & OLS Regressions  

###All villages (dat1ineligible)  ##No fixed
alliv<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+ hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+ hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=dat1ineligible)
summary(alliv)
se.func(alliv, dat1ineligible)

nrow(dat1ineligible)  ##N is 4,211
mean(dat1ineligible$sc)  # dep var mean is .587



####### Looking at the villages below .67 of fraction eligible  (nothing significant here.  Many things zero)
low.5OLS.A <- lm(sc~Trpool, data=lowelig.5)
summary(low.5OLS.A)
se.func(low.5OLS.A, lowelig.5) 


lowelig.OLS.B <- lm(sc ~ Trpool + gender + indig + welfare + hohedu + hohage + hohsex + fam_n+dist_sec, data=lowelig.5)
summary(lowelig.OLS.B)
se.func(lowelig.OLS.B, lowelig.5)


nrow(lowelig.5)
mean(lowelig.5$sc)

#w/ state fixed effects and village fixed
lowelig.OLS.B2 <- lm(sc ~ Trpool + gender + indig + welfare + hohedu + hohage + hohsex + fam_n + dist_sec +Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=lowelig.5)
summary(lowelig.OLS.B2)
se.func(lowelig.OLS.B2, lowelig.5)



###Subsetting for 98  
loweligible98 <- lowelig.5[lowelig.5$year==98,]
head(loweligible98)


lowelig.OLS.98<- lm(sc ~ Trpool + gender + indig + welfare + hohedu + hohage + hohsex + fam_n + dist_sec, data=loweligible98)
summary(lowelig.OLS.98)
se.func(lowelig.OLS.98, loweligible98)



###w/state fixed effects
lowelig.OLS.98.2<- lm(sc ~ Trpool + gender + indig + welfare + hohedu + hohage + hohsex + fam_n + dist_sec + factor(state)-1, data=loweligible98)
summary(lowelig.OLS.98.2)
se.func(lowelig.OLS.98.2, loweligible98)




   
###Subsetting for 99
loweligible99 <- lowelig.5[lowelig.5$year==99,]
head(loweligible99)


lowelig.OLS.99<- lm(sc ~ Trpool + gender + indig + welfare + hohedu + hohage + hohsex + fam_n+dist_sec, data=loweligible99)
summary(lowelig.OLS.99)
se.func(lowelig.OLS.99, loweligible99)



###w/state fixed effects
lowelig.OLS.99.2<- lm(sc ~ Trpool + gender + indig + welfare + hohedu + hohage + hohsex + fam_n + dist_sec + factor(state)-1, data=loweligible99)
summary(lowelig.OLS.99.2)
se.func(lowelig.OLS.99.2, loweligible99)








##IV ###Estimates are usually positive, but relatively small (smaller than what they find), insignificant, and sometimes negative -> AMBIVALENT
###w/no fixed effects
table5A4.lowelig<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+ hohsex+fam_n+dist_sec|Trpool+gender+indig+welfare+hohedu+hohage+ hohsex+fam_n+dist_sec, data=lowelig.5)
summary(table5A4.lowelig)
se.func(table5A4.lowelig,lowelig.5)               



### w/ fixed effects by state 
table5A4.lowelig2<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+ hohsex+fam_n+dist_sec+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+factor(state)-1, data=lowelig.5)
summary(table5A4.lowelig2)
se.func(table5A4.lowelig2,lowelig.5)



###W/ state fixed effects and village controls 
table5A4.lowelig3<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+ hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=lowelig.5)
summary(table5A4.lowelig3) 
se.func(table5A4.lowelig3,lowelig.5)


###FIRST STAGE
lm1 <- lm(Esc10r~Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=lowelig.5)
se.func(lm1, lowelig.5)



#W/ only village level effects (no state)
table5A4.lowelig4<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+ hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r, data=lowelig.5)
summary(table5A4.lowelig4)
se.func(table5A4.lowelig4,lowelig.5)





#####IV by year  ##98 not significant and small
table5A4.lowelig98 <-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=loweligible98)
summary(table5A4.lowelig98)
se.func(table5A4.lowelig98,loweligible98)

##FIRST STAGE
lm2 <- lm(Esc10r~Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=loweligible98)
se.func(lm2, loweligible98)


### Year 99 ## also negative and insignificant
table5A4.lowelig99 <-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r + factor(state)-1, data=loweligible99)
summary(table5A4.lowelig99)
se.func(table5A4.lowelig99,loweligible99)


##FIRST STAGE
lm3 <- lm(Esc10r~Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r + factor(state)-1, data=loweligible99)
se.func(lm3, loweligible99)





###By grade  ###Significant, STRONG effect for just 6th and 7th
###NOTE- might want to check with and without fixed effects to see if effect is still there
lowelig.67 <- ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+gender+indig+year+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+
Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=lowelig.5[lowelig.5$grade97<7,])
summary(lowelig.67)
se.func(lowelig.67, lowelig.5[lowelig.5$grade97<7,])


###The coefficient is negative and insignificant for 8th and 9th
lowelig.89 <- ivreg(sc~Esc10r+gender+indig+year+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+year+gender+indig+year+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=lowelig.5[lowelig.5$grade97>6,])
summary(lowelig.89)
se.func(lowelig.89, lowelig.5[lowelig.5$grade97>6,])


###Looking at just the kids below the median welfare (of 900)
###Large, but insignificant effect
lowelig.lowwelfare<-ivreg(sc~Esc10r+year+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|
	Trpool+year+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=lowelig.5[lowelig.5$welfare<=900,])
summary(lowelig.lowwelfare)
se.func(lowelig.lowwelfare,lowelig.5[lowelig.5$welfare<=900,])








##############################################################
###### Looking at villages ABOVE .6667 fraction eligible.  

highelig.OLS.A <- lm(sc~Trpool+ gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=highelig.5 )  
summary(highelig.OLS.A)
se.func(highelig.OLS.A, highelig.5 )


nrow(highelig.5)
mean(highelig.5$sc)

highelig.OLS.B <- lm(sc ~ Trpool + gender + indig + welfare + hohedu + hohage + hohsex + fam_n + dist_sec, data=highelig.5)
se.func(highelig.OLS.B, highelig.5)


 

###Subsetting for 98
higheligible98 <- highelig.5[highelig.5$year==98,]
head(higheligible98)


highelig.OLS.98<- lm(sc ~ Trpool + gender + indig + welfare + hohedu + hohage + hohsex + fam_n + dist_sec, data=higheligible98)
summary(highelig.OLS.98)
se.func(highelig.OLS.98, higheligible98)


###Subsetting for 99
higheligible99 <- highelig.5[highelig.5$year==99,]
head(higheligible99)



#no state effects
highelig.OLS.99.1<- lm(sc ~ Trpool + gender + indig + welfare + hohedu + hohage + hohsex + fam_n + dist_sec, data=higheligible99)
se.func(highelig.OLS.99.1, higheligible99)


#W/ state effects
highelig.OLS.99<- lm(sc ~ Trpool + gender + indig + welfare + hohedu + hohage + hohsex + fam_n + dist_sec+ factor(state)-1, data=higheligible99)
se.func(highelig.OLS.99, higheligible99)

nrow(lowelig.5) # N is 2153
mean(lowelig.5$sc) #Dep var mean is .584
nrow(loweligible98) # N is 1081
mean(loweligible98$sc) #Dep var mean is .654
nrow(loweligible99) # N is 1072
mean(loweligible99$sc) #Dep var mean is .513


###IV ###  
## Village and state fixed effects  ##LARGER EFFECT THAN THEY FOUND.  HIGHLY SIGNIFICANT
table5A4.highelig<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=highelig.5)
summary(table5A4.highelig)
se.func(table5A4.highelig,highelig.5)


#FIRST STAGE
lm4 <- lm(Esc10r~Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=highelig.5)
summary(lm4)
se.func(lm4, highelig.5)


### No fixed effects
table5A4.highelig2<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec, data=highelig.5)
summary(table5A4.highelig2)
se.func(table5A4.highelig2,highelig.5)




### Only village fixed effects
table5A4.highelig3<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r, data=highelig.5)
summary(table5A4.highelig3)
se.func(table5A4.highelig3,highelig.5)



### Only state fixed effects
table5A4.highelig4<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+factor(state)-1, data=highelig.5)
summary(table5A4.highelig3)
se.func(table5A4.highelig4,highelig.5)






#####IV by year  ###BIGGER effect for 98
table5A4.highelig98<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=higheligible98)
summary(table5A4.highelig98)
se.func(table5A4.highelig98,higheligible98)

##FIRST STAGE
lm5 <- lm(Esc10r ~ Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=higheligible98)
se.func(lm5, higheligible98)


### Year 99
table5A4.highelig99<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=higheligible99)
summary(table5A4.highelig99)
se.func(table5A4.highelig99,higheligible99)


#FIRST STAGE
lm6 <- lm(Esc10r~Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=higheligible99)
se.func(lm6,higheligible99)


nrow(highelig.5) ##N is 2,058
mean(highelig.5$sc) # Mean of dep var is .591
nrow(higheligible98) #N is 1,035
mean(higheligible98$sc) # Mean of dep var is .662
nrow(higheligible99) #N is 1,023
mean(higheligible99$sc) # Mean of dep var is .519

###By grade 6th and 7th ###positive, insignificant
highelig.67 <- ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=highelig.5[highelig.5$grade97<7,])
se.func(highelig.67, highelig.5[highelig.5$grade97<7,])


length(unique(dat1ineligible$id)) ##2446 kids ineligible (NOT OBS)
nrow(dat1) ##21,705 obs total
nrow(dat1[dat1$eligible=="pobre",]) ###17494 eligible obs
length(unique(dat1$id[dat1$eligible == "pobre"]))  ##10024 eligible kids (NOT OBS)

###8th and 9th
##### Sign turns negative, insignificant
highelig.89 <- ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=highelig.5[highelig.5$grade97>6,])
se.func(highelig.89, highelig.5[highelig.5$grade97>6,])





###Looking at just the kids below the median welfare (of 900)
###sizable effect that are significant
highelig.lowwelfare<-ivreg(sc~Esc10r+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1|Trpool+gender+indig+welfare+hohedu+hohage+hohsex+fam_n+dist_sec+Egender10r+Eindig10r+Efam_n10r+Ehohedu10r+Ehohage10r+Ehohgender10r+factor(state)-1, data=highelig.5[highelig.5$welfare<=900,])
summary(highelig.lowwelfare)
se.func(highelig.lowwelfare,highelig.5[highelig.5$welfare<=900,])



